A model for the transport of muons in extensive air showers 



L. Cazon a *, R. Conceicao a , M. Pimenta a ' b , E. Santos a 

"LIP, Av. Elias Garcia, 14-1, 1000-149 Lisboa, Portugal 
b Departamento de Fisica, 1ST, Av. Rovisco Pais, 1049-001 Lisboa, Portugal 



Abstract 

In this article we identify the key elements that govern the propagation of muons from the produc- 
tion in extensive air showers to ground. We describe a model based on simple assumptions that 
propagates the muons starting from the few relevant distributions at production. We compare the 
results to the ground distributions given by a full air shower Monte Carlo. This study is motivated 
by the need of modeling the muon component in extensive air showers with the goal of experi- 
mentally reconstructing their distributions at production, which act as a footprint of the hadronic 
cascade. 
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1. Introduction 

The nature of the Ultra High Energy Cosmic Rays remains unknown. The state of the art ex- 
periments have not yet understood key aspects necessary to answer this question. While the energy 
and arrival direction of the cosmic rays to Earth can be fairly well reconstructed, the primary mass 
is difficult to determine. The high energy spectrum may be fitted by a number of combinations of 
light or heavy nuclei, since the density evolution and maximum energy achievable at the sources 
are yet unkown. The observation of anisotropics on the cosmic rays sky does not necessarily favor 
either light or heavy nuclei, because the galactic and extragalactic magnetic fields and the location 
of the sources are still unknown. The direct mass determination from air shower observables is 
not conclusive either. Our understanding of the hadronic particle physics is only supported up to 
the LHC energies. Beyond those energies we must rely on the extrapolations of the hadronic in- 
teractions models which diverge from one to another. Besides, these new kinematic regions might 
uncover new phenomena not yet accommodated in models. Changes in the hadronic physics and 
in the composition of the primary share a region of the phase space, being difficult to break the 
degeneracy and answer both questions. 

Recent results from the Pierre Auger Observatory JH] on the evolution of the depth of the elec- 
tromagnetic shower maximum have been usually interpreted as a change towards heavier composi- 
tion at the highest energies, provided that the extrapolations of the hadronic interaction models are 
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correct. Nevertheless, an abrupt change on the hadronic interactions at the highest energies could 
be possible, leading to a rapid increase of the cross section yj] and changes on other aspects of the 
multiparticle hadronic production. Moreover, recent results from the LHC indicate that the current 
understanding of the forward direction embedded on the hadronic interactions models might be 
insufficient flU. 

Auger has also shown [0] that the number of muons in Extensive Air Showers (EAS) is under- 
estimated by the current hadronic interactions models, even for the case of iron primaries by 40% 
with respect to QGSJET-II [0,0]. Muons in extensive air showers have been already the subject 
of many experimental studies, from their absolute number at ground to the longitudinal profile 
at different energies. A book by Grieder [8I] contains an excellent compilation of most available 
results. Very recently, the KASCADE-Grande collaboration has reported the measurement of the 
longitudinal profile of the production points of muons by tracking the trajectories of the detected 
muons at ground back to the shower axis DSD. They have also published the number of muons at 
ground for showers with different electron richness [10]. KASCADE-Grande reaches up to 10 18 
eV. 

Despite of being a detector that was not originally optimized for muon reconstruction, the 
Auger Collaboration has recently measured 111 ill the maximum of the muon production depth 
profiles at energies above 10 19 2 eV by mapping the arrival time of muons far from the core onto 
muon production distances ill 211 . 

Recent plans from Auger include the deployment of a set of buried detectors aimed to the 
muonic component [13], allowing us to explore with low systematics the region below 10 18 eV, 
which overlaps with KASCADE-Grande and also with the LHC results. 

This work is motivated by the need to extract complementary information carried by the 
muons, truly hadronic messengers in EAS. We aim to properly model the mechanisms that govern 
the muon distributions in air showers in order to peer into the details of the hadronic shower. 

EAS develop in a complex way as a hadronic multiparticle production that generates a hadronic 
and an electromagnetic cascade which travel down the atmosphere. The electromagnetic (EM) 
component spreads out in time and space, reaching more than 1 km away from the shower core in 
large numbers, enough for detectors placed at ground level (several square meters of surface) to 
record their signal. On the other hand, given their large density near the core, they produce sizable 
amounts of fluorescence light that can be detected with far away ultraviolet telescopes, recording 
their longitudinal development in moonless nights. 

The hadronic cascade is much less numerous and thus hard to detect directly. It consists mainly 
of low energy pions, and fewer high energetic particles, such as leading baryons and mesons that 
carry a large fraction of the primary energy deep into the atmosphere. The energy and momentum 
of these leading particles depend on the details of the high energy hadronic interaction models, 
which determine the production of the lower energy bulk of mesons through the inelasticity and 
multiplicity of the interactions. The hadronic cascade is the main engine of the air shower: it feeds 
the EM cascade mainly by the decay of neutral pions and also feeds back through the interaction 
of charged pions, or by means of the less numerous kaons. When a pion or kaon decays into 
a muon, the muon might leave the hadronic core and transports information far away from the 
central reg ion. 

In Ill4lll5ll . and later updated in II 1611 . it was shown that the arrival time distributions of muons 
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at ground emerge as a direct transformation of the muon production depth distribution and the 
energy spectrum at ground. In the present paper we develop the model that links all the relevant 
distributions at the moment of production, namely, the muon energy, transverse momentum and 
production depth to the observed muon distributions at ground, namely the energy spectrum, the 
arrival time delay distribution, the apparent production depth distribution and the lateral distribu- 
tion (that is, the muon surface density at ground) Q. This knowledge is useful for fast air shower 
Monte Carlo simulations like CONEX y/70, and for use in the reconstruction algorithms of the 
number of muons, muon production depth, production energy and possibly the transverse momen- 
tum distributions, all of them directly inherited from the hadronic cascade. The distributions of 
all muons at the moment of production, that is, including those that would decay later on flight, 
exhibits the most universal features independently of the observational conditions. 

The inference of the fundamental distributions of all muons at production from the observed 
distributions it is not straightforward, since the information carried by muons below some energy 
threshold - defined by the amount of matter transversed in the atmosphere - is completely lost. 
Nevertheless, by combining different observation conditions (zenith angle and distance to core) 
which correspond to different effective energy thresholds, one might constrain the average distri- 
butions at production. Notice also that the propagation of muons itself does not depend on the 
details of the hadronic interactions models, being a completely decoupled problem. 

A detailed understanding of the detector response to the different particles of an air shower is 
necessary to properly identify and interpret the muon information. Such a study is out of the scope 
of this paper. 

This paper is organized as follows. In section 2 we define some quantities and describe the 
muon distributions at production. In section 3 we describe the propagation of muons and the 
approximations used. In section 4 we analyze the distributions at ground after propagation and 
compare them to a full air shower Monte Carlo simulation. In section 5 we discuss the effects 
of averaging the energy and transverse momentum distributions. Finally, we comment on the 
prospects and conclusions. 

2. The production of muons in EAS 

When a cosmic ray enters the atmosphere it creates an air shower of particles. The extrapola- 
tion of the original trajectory defines the so called shower axis and its intersection with the ground 
surface defines the shower core. We will use a Cartesian coordinate system which is centered at 
the core position in ground, with the z-axis parallel to the shower axis. The y axis will be parallel to 
ground, and the x-axis is positive downwards, entering the earth with an angle 9 with the surface, 
which is also the angle between the z-axis and the zenith's direction (see Fig. [TJ)- A cylindrical 
coordinate system can be defined by r = ^Jx 2 + y 2 and £ = arctanj/x and it will be sometimes 
used for convenience. 



Another observable at ground is the angular distribution, which allows the reconstruction of the apparent produc- 
tion depth distribution by backtracking the muon trajectories to the shower axis. This will be studied in detail and will 
be published elsewhere B26I1 . 
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Figure 1 : Scheme showing the shower plane (perpendicular to the shower axis), the ground surface and the system of 
coordinates. 



The z-coordinate can be expressed in terms of the amount of matter along the shower axis from 
the top of the atmosphere: 

P(z')dz'. (1) 

We will use indistinctly z and X to express the position along the shower axis where muons are pro- 
duced. The z variable is more suited for calculations regarding geometry or kinematics, whereas 
the X variable is used for the evolution of the cascade. 

In Ill4ll . it was argued that the transverse position of the production of muons, thus of the 
parent mesons decay, is confined to a relatively narrow cylinder. Fig. [2] right panel, displays 
the distribution of the y-coordinate where muons were produced (the shower axis is at y=0) for 
different primaries. On the right panel, the average and the v-coordinate containing 50% and 
90% of the production points are displayed as a function of the atmospheric depth. The average 
value is of tens of meters. This distance is small when compared to the distances involved in 
EAS experiments, which span from hundreds of meters to several kilometers in the perpendicular 
plane. For instance, the Pierre Auger Observatory has its tanks separated by 1 .5 km. Therefore, 
the position where the muon has been produced can be approximated by (0, 0, z), or simply z. 

Every dX along the shower axis, dN muons are produced within a given energy and transverse 
momentum interval dE t and dp,. Their overall distribution at production can be described in 
general with a 3-dimensional function, as: 

d 3 N 

————= F(X,E h cp t ) (2) 
dX dhjdcpt 

These 3D-distributions were recorded during simulation with CORSIKA v6.980 lfl8ll at the mo- 
ment of production, along with the standard ground particle output files. These distributions define 
most of the knowledge about muons contained in the air showers. 
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Figure 2: Left Panel: average distribution of the positive y-coordinate at the production point of muons for different 
primaries at 10 19 eV at 60 deg. Right panel: average, median and 90% quantiles of the y distribution for different 
depths for a proton shower of 10 19 eV at 60 deg. 
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Figure 3: Left panel:total number of muons produced per g cirT 2 (h(X)) for 50 proton showers at 10 19 eV and 60 deg. 
Right Panel: normalized spectrum of muons at production at two distances to the shower maximum, X' - X - X^, ax 
300 g cirT 2 and X' - X - X^. dx = 300 g cirT 2 , for the same showers. 



A library of CORSIKA v6.980 showers was created, with samples of 50 proton showers at 0, 
40, 60 and 70 degrees with QGSJET-II.03 [0,0] model and energy 10 19 eV. The relative thinning 
was set to 10" 6 , the maximum weight was set to 10 4 and the hadronic maximum weight of 10 2 . 
The inner radial thinning was set to 1 cm. The kinetic energy cuts where set to 0.05 GeV for 
muons and hadrons and 0.003 GeV for electrons and photons. The altitude was set to 1400 m a.s.l, 
with the magnetic field of Malargue, Argentina. 
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At 60 degrees, we also run 50 shower samples of iron and photon primaries with QGSJET- 
11.03, and for iron and proton for the models SIBYLL2.1 0] and EPOS 1.99 HQ. We also 
run proton showers at 10 18 5 eV, 10 19 eV and 10 19 5 eV with QGSJET-II.03. 

A sub-sample of proton QGSJET-II.03 10 19 eV showers was simulated with a lower maximum 
weight of 10 3 for EM particles and 10 for hadrons (and muons). This subset was used to check the 
effects of the thinning, which are negligible. 

The projection into the X (or z) axis becomes 



h(X) = J F(X, E h cp t )dEidcp t (3) 

and it is the so called total/true Muon Production Depth (Distance) distribution, or MPD-distribution 
for short. It does not depend on the observational conditions since it does not contain any prop- 
agation effects of muons through the atmosphere. Notice that this is different from the MPD- 
distributions of detected muons at a given position on ground 4f which includes the effects of 
propagation, as it will be explained later. This distribution is sometimes referred to as apparent 
MPD-distribution. 

The total number of muons produced in a shower is 



N = j h(X)dX (4) 

It should be noted that this number is intrinsically different from the surviving muons, which 
are affected by the fluctuations of the depth of the first interaction, which changes the distance 
traveled by muons to ground. Some of the techniques used by Auger I122ll use a fixed distance to 
the shower core, so they can also be affected by the lateral spread of the parent mesons. Notice that 
in CORSIKA, the function F(X, E h cp t ) is only known above a certain energy threshold, E t > E th . 
Therefore the total/true MPD-distribution depends on E th , becoming h Eth (X). In the same manner, 
the total number of produced muons corresponds to the total number of muons above the energy 
threshold, Ne,,,- The simulations of the present work have a total energy threshold E,/, = 0.155 
GeV, the lowest allowed by CORSIKA code. The differences in the distributions at observation 
level induced by this particular threshold are negligible because such low energy muons decay on 
flight before reaching ground. E th value must be specified when referring to Ne, h or h Eth for which 
the low energy muons can contribute in large amounts B23I1 . Changing from a concrete value of E th 
into another is trivial provided that we know F(X, Ej, cp t ). 
Eq. [2] can be factorized and expressed as the product 

F{X, E u cp t ) = h(X) f x (Ei, cp t ) (5) 

where the function f x (Ei,cp t ) = F(Z ^' ) cp ' ) becomes the normalized Ej and cp, distribution at a 
given production depth X. 



In the approximations made on Hi 41 LL5_L LLmL fx did not depend on X and it was factorized in 
2 independent distributions on E, and cp t . This allowed analytical approximations of the distribu- 
tions at ground. In this work we have included these correlations, improving the accuracy of the 
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energy, production depth, and time distributions at ground, and allowing for a proper description 
of the muon lateral distribution at ground. 

The function h(X) tracks the longitudinal development of the hadronic cascade and represents 
the production rate of muons per gem" 2 . Its shape and features are extensively discussed in ll23ll . 
The depth at which h(X) reaches a maximum is denoted as X^ ax . X^ ax correlates with the first 
interaction point Xi which corresponds to the first interaction of the primary in the atmosphere and 
the start of the cascading process I123n . The most important source of fluctuations in air showers 
corresponds to the fluctuations of X\, which causes an overall displacement of the whole cascade 
at first approximation. The amount X' = X - X^ ax defines the amount of matter with respect to the 
shower maximum. The distributions can be expressed in terms of X', where the most important 
source of fluctuations has been eliminated, and only the remaining effects are present. 

In Fig. [3] (left panel) h(X) is shown for a sample of 50 showers. The fluctuations on the 
normalization and on X^ are clearly observed. In the right panel we can see the normalized 
energy spectrum for two values of X', namely -300 g cm" 2 and 300 g cm" 2 . Both the energy and 
the transverse momentum show similar features when referred to the same distance to the shower 
maximum, X' . From now on, whenever we average distributions at production we do it on X' , that 
is, matching the maxima. 

In IU4I1 and lllal the muon spectrum at production was approximated by a power law, Et 2 - 6 , 
following the high energy tails of the pion production on the hadronic reactions. In that approach 
we did not have access to all muons at production, and the energy spectrum was extrapolated down 
to low energies with the same power law. In Fig. H left panel, the actual average energy spectrum 
of all muons at production is displayed for proton showers at 10 19 eV in different X' layers. At 
low energies the single power law clearly does not work. In addition, the energy spectrum evolves 
with X' by becoming softer, and stabilizing the shape after the shower maximum. In Fig. |5l 
left panel, the energy spectrum is displayed for different zenith angles at X' = 0, showing a mild 
dependence, becoming harder at higher zenith angles. This might be due to the higher critical 
energy of the pions at higher zenith angles, given that the shower develops in less dense air in 
average. Regarding the dependence with the primaries, Fig. [51 right panel displays the spectrum 
for proton, iron and photon primaries at X' = 0. Whereas proton and iron curves practically 
overlap, photons, on the other hand, produce a much softer energy spectrum. The mechanisms for 
production of muons in photon showers is basically photo-pion production, so they are intimately 
related to the energy spectrum of the EM cascade. Lastly, Fig. [6j left panel, displays a comparison 
between different models, in the same conditions. While QGSJET- 11.03 and SIBYLL2.1 overlap, 
EPOS 1 .99 shows a slightly different behavior, with a softer spectrum. Fig. [6l right panel, displays 
a comparison of different primary energies with the same hadronic interaction model, where the 
curves practically overlap. 

The transverse momentum distributions are responsible for most of the lateral displacement 
of muons with respect to the shower axis. In [15], the p t distributions were approximated by 
an unique function, dN/dp t = Pt/Q 2 exp(-p t /Q), independent of the energy of the muon and 
its production depth, primary mass and zenith angle. In the current work, we uncover in detail 
all the dependencies. The p, distributions display a quite universal shape as a function of the 
zenith angle, (see left panel of Fig. [7] displaying the normalized average distributions of the p t 
distributions at X' = g cm" 2 ), and do not depend on the primary (right panel Fig. [7]), except 
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Figure 4: Normalized average energy distribution of all muons at production for proton showers at 10 19 eV and 60 deg 
zenith angle simulated with QGSJET-II.03. Left panel shows cuts at different X' layers. Right panel: 50% quantile 
(median) as a function of X' along with a band defined by the 16% and 84% quantiles. The average value is also 
plotted as a solid line. 



9 = deg 




— proton 




log l0 (E/GeV) 



Figure 5: Normalized average energy distribution of all muons at production for proton initiated showers at 10 eV 
simulated with QGSJET-II.03for different zenith angles (left panel), and different primaries at 60 deg zenith angle 
(right panel). 



for the case of photons, which respond to a completely different muon production mechanism, as 
explained before. As the shower evolves, the p t spectrum becomes softer (Fig. [81 left panel shows 
the evolution as a function of X'). Besides this dependence on X', the p, distributions also depend 
on the energy of the muons. The right panel shows the medians of the ^distribution for 2 different 
energy cuts. The low energy muons display a smaller p t , and at high energies, the p t distribution 
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Figure 6: Normalized average energy distribution of all muons at production for proton showers at X' = g cirT 2 and 
60 deg zenith angle. Left panel, 10 19 eV for different hadronic interactions models. Right panel, QGSJET-II.03 for 
different primary energies. 



prefers higher p t values. Both distributions show a different dependence on X'. We have found 
that the different correlations of the p t with E t and X must be included into the model in order to 
properly predict the muon lateral distribution at ground. 

Fig. |9l left panel, displays the p, distribution for 3 models. The high p t tail of SIBYLL2.1 is 
suppressed with respect to the other models, which practically overlap. Fig. [9} right panel, displays 
a comparison of different primary energies with the same hadronic interaction model, where the 
curves practically overlap. 



3. Propagation of muons through the atmosphere 



The propagation of muons in matter below 1 TeV is reviewed in [25]. In this section we 
will describe the aspects relevant for the distributions of muons at ground in UHECR-induced air 
showers. The extremely energetic muons coming from the first interactions would need a huge 
area covered by high resolution detectors in order to collect and identify them. On the contrary, 
the region below 1 TeV constitutes the bulk of the muons, and it is the most important from the 
statistics point of view. 

The transport of muons to ground was implemented by means of a fast Monte Carlo in two 
steps. In the first step the muons were propagated to ground following a straight line according to 
their 3-momentum and the continuous energy loss was calculated. In the second step, the multiple 
scattering and magnetic field effects were included, the impact point on ground was corrected, and 
the energy loss reevaluated. Then, the probability of decay and the time delay for the corrected 
energy loss and trajectory were obtained. 

For each shower, the variables at the production point, namely z, E h and cp t were sampled 
from the distributions F(z, E t , cp t ) with a weight w, . The momentum of the muon defines an angle 
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Figure 7: Normalized average p, distribution of all muons at production for showers at 10 eV at X' — gem - 
simulated with QGSJET-II.03. Left panel, different zenith angles, with proton primaries. Right panel, different 
primaries at 60 deg zenith angle. 
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Figure 8: Normalized average p, distribution of all muons at production for proton initiated showers at 10 eV 
simulated with QGSJET-II.03 for different values of X' (left panel). Median cp, as a function of X' for the total 
number of muons, and the median for those muons with Ej > 2 GeV and E, < 2 GeV, as labeled (right panel). 



a respect to the shower axis given by 

cp t 

sina^— , (6) 

Ei 

where we took the approximation cp - E t . A 10 GeV (1 GeV) muon typically will span a 1 
deg (10 deg) outgoing angle with respect to the shower axis. The outgoing polar angle follows a 
symmetric distribution ^ = ^ that was also randomly sampled. 
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Figure 9: Normalized average p, distribution of all muons at production for proton showers at X' = gcm~ 2 and 
60 deg zenith angle. Left panel, 10 19 eV for different hadronic interactions models. Right panel, QGSJET-II.03 for 
different primary energies. 



From the production point (0, 0, z) a straight line trajectory defined by (a, £) is extrapolated 
until it hits ground, defining the arrival position of the muon (r, g, A), where A is the z coordinate 
of the ground in the shower system, and r is the distance to the core in the perpendicular plane. 
For a flat ground surface they become 

r = t- (7) 

cos £ tan 9 + — — 

' tana 

A = r cos 4" tan (8) 

which can be generalized for a curved earth surface. The distance traveled by the muon from the 
production point to the ground is 

I = V>" 2 + (z - A) 2 (9) 



As stated in 111411 . we can define a plane parallel to the xy plane that travels at the speed of light 
that contained the first interaction point and hits ground at t = 0. Assuming that all particles travel 
at speed of light and no delay is accumulated during the hadronic cascade, the muon is produced 
with no delay. The arrival time delay of the muon respect to this plane front becomes 

ct g = l-(z-A) (10) 

which can be approximated by ct g - in most practical cases. A correction due to the path 
traveled by the parent meson was applied. The trajectory would have started an amount Az n = 
ct^^j cos a higher up in the atmosphere. Thus, we replace z (and /) in Eq. Q/Jby z + Az n . A muon 
produced at z = 10 km from ground and reaching r = 1000 m will be delayed t g - 165 ns because 
of geometric effects, (t g - 167 ns if we did not take the parent meson correction). 
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Along this trajectory the muon will suffer a number of processes, namely: energy loss, multiple 
scattering, and magnetic field deflections. Each one of these processes will modify the momentum, 
trajectory and time delay with respect to the plane front traveling at velocity c. In addition, muons 
can decay on flight with a certain probability, in which case some information is lost. 

3.1. Continuous energy loss in the atmosphere 



Muons lose energy as they travel through the atmosphere [25]. Only the energy loss by ioniza- 
tion is taken into account as: 

dE 

-a (11) 



and where a has been parametrized as 

a = 1 2.06 + 0.5453^ + 



dX 



0.0324 



x 10" 3 (GeV/gcnT z ) 



-2\ 



(12) 



(£ + 1.0312) 2 , 

where £ = log 10 (£7GeV). For energies E >~ 50 GeV, where the radiative losses start to appear, 
the relative error on the energy loss to the total energy are neglected. 

X(z) can be calculated with the detailed description of the atmosphere in different layers cur- 
rently used by most Monte Carlo codes. Nevertheless, we used a compact description of the 
atmosphere in order to calculate the energy losses of muons in a single step. We used the expo- 

_ h 

nential approximation in a single layer, p = p e h » . The thickness of atmosphere traversed by a 
muon between the production point and ground becomes 



AX(l) = p(l)dl = pok [l-e 



i) 



(13) 



where / = h/ cos 8^ and l Q = ho/ cos 8^. h is the height and and 8^ is the zenith angle of the 
trajectory of the muon. The values of po and h were obtained by the best fit to X(z) up to 10 km 
height. Notice that the region close to ground is the one that affects the most the propagation of 
muons. 

The energy of the muon evolves as 



E(l) = Ei-apok^ -e '») 



(14) 



The finite energy of the muons induces a delay with respect to a particle traveling at the speed 
of light. The particles belonging to the hadronic shower decrease their energy on a geometric 
progression, due to the cascading process. The delay accumulated prior to the moment of decay 
into a muon is neglected. After the muon is produced, its kinematic delay becomes: 



ct f 



1 - 



mc 

W) 



dl' -I 



(15) 



This equation can be easily integrated for the case where ^ < 1 yielding: 



ImV 
ct e = -— — 



£(/) 

,1 1 \ I Z ologf 
-la —-—+ — + — = — - 



E, 



(16) 
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where 

E f = E i -ap Q l Q (l-e-^ (17) 

and 

E co = E f +p al (18) 

Most of the delay of the particle is accumulated at the end of the trajectory, where the energy is 
lower due to the continuous energy loss. A 5 GeV (10 GeV) muon traveling 10 km with a zenith 
angle of 60 degrees, will reach ground with an energy Ef = 3.0 GeV (7.8 GeV) and a kinematic 
delay of 12 ns (2.3 ns). 

3.2. Probability of Decay 

Muons can decay on flight with a probability that depends on the energy. The low energy 
muons are reduced or totally suppressed if they are below a certain energy threshold dependent on 
the amount of matter they would have to traverse to ground, which is of the order of 0.2 GeV every 
~ 100 gem -2 . The probability of decay in a interval dl is -—j?dl. By integrating, we obtain the 
probability of decay as a function of the traveled distance: 

mc 2 '0 

_L1 g — KT (19) 

In our Monte Carlo, the weight of the muon at ground is substituted accordingly 

w = w i p(E i ,l) (20) 

A 5 GeV (10 GeV) muon traveling 10 km to ground with a zenith angle of 60 degrees has a 
probability of survival p = 0.67% (0.84%). 

Decay plays a fundamental role on shaping the distributions of muons at ground: the muon 
lateral distribution, time distribution and the energy spectrum are suppressed in the regions domi- 
nated by low energy muons, if we compared to the case where the decay was not present. 

3.3. Geomagnetic Field 

The magnetic field of the Earth affects the muons by bending their trajectory, changing the 
impact point on ground and delaying the arrival time. The effects are more visible the longer the 
trajectories. Muons travel long paths in inclined showers, reaching up to 220 km for showers at 
86° zenith angle. In 1E4I1 a model that describes this effect was developed. After neglecting the 
effects of the transverse momentum, cp t <£we obtain the radius of curvature R of the trajectory 
of a muon of energy E : 

R = J^^^- (21) 
eB ± ceB ± 

where B ± is the projection of the magnetic field onto the perpendicular plane. The energy of the 
muon was approximated by E = {Ei + E/)/2. 
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Figure 10: Deviation from the expected trajectory of a positive particle (traveling downwards) because of the magnetic 
field (entering the paper). The difference between the arc s and a chord /' introduces a time delay. 



Let us now define a set of coordinates (x B , y B ) in the shower plane in such a way that the positive 
y B axis is aligned with B ± . As a muon travels it will be shifted 8x B in the direction perpendicular 
to B ± given by (see Fig. [TT 



6x R =R 




1Z 2 

-2R (22) 



The + sign depends on the charge of the muon, which is sampled randomly with equal probability. 
This means that the locations of arriving muons are shifted in dependence on the traveled distance 
/ but also on the energy of the muon. The coordinate y B is left unchanged. For a typical value of 
B ± = 10 f/T, a 5 GeV (10 GeV) muon that travels 10 km deviates 5x B = 83 m (16.7 m). 
The actual path length of the trajectory is enlarged as 

s(l) = l + 2R arcsin ^ =* I + (23) 

The difference between s(l) and / is the so called geomagnetic time delay: 

1 I 3 

A 5 GeV (10 GeV) muon traveling 10 km, will delay t B 0.04 ns (0.01 ns). 

3.4. Multiple Scattering 

Muons traversing the air are deflected by small-angle scatters. Most of these deflections are 
due to Coulomb scattering from nuclei. The many small interaction add up to an angular and space 
deviation that follows a Gaussian distribution to a good approximation. 
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Figure 1 1 : Examples of RMS displacement due to multiple scattering for different zenith angles and Ef = 1 Ge V and 
Ef = 10 GeV as labeled. The solid black line is the parametrization given in l25ll with E = (£, + £/)/2 



Using a dedicated library of CORSIKA simulations with the geomagnetic field off, we have 
created look-up tables (in the form of collections of 3D histogram) with the different relevant 
scattering quantities, namely, the RMS of the lateral displacement, the RMS and mean of the 
logarithm of the time delays, and the RMS and mean scattering angle as a function of the energy 
at ground and total distance traveled by the muon. The lateral displacement was calculated by 
comparing the impact point the muon would have if we extrapolated its momentum at production 
by a straight line trajectory with the actual impact point at ground. The time delay was calculated 
by subtracting the known geometric and kinematic delays from the total delay given by CORSIKA. 
This remaining time delay includes the effects of the multiple scattering and also the differences 
between the real geometric and kinematic delays and the approximations made in the previous 
sections. The scattering angles were calculated by comparing the momentum at production with 



the momentum at ground, and will be used in [26]. 



Fig. QT] displays the RMS of the lateral displacement as a function of the traveled depth for 
different zenith angles with a final energy 1 and 10 GeV. The solid black line represents the analytic 



approximation given in 02511 . The effect of the multiple scattering on the position of the muons is 
accounted for in the model Monte Carlo as follows: the position of the muon is displaced from 
its geometric extrapolation following a 2-Gaussian distribution with RMS given by the look-up 
tables. The arrival time delay is also corrected by adding the extra term that accounts for the 
multiple scattering and the precision effects, sampled from a log-Gaussian whose parameters were 
also recorded in the look-up tables. 

After the geomagnetic and multiple scattering correction, the total distance I from the produc- 
tion point is reevaluated and the final energy is recalculated. The effects of the multiple scattering 
on all the analyzed profiles, where the distances to the core are larger than 100 m, were found to 
be negligible. This is so because in most practical cases the width of the profile is wider than the 
width introduced by the multiple scattering smearing. Notice that when analyzing muon by muon 
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these effects might become important in some cases. 



4. Ground distributions 

The model Monte Carlo transports A/b/w* muons from the production point to the ground level. 
The weight at production w, can be selected as a trade off between speed and statistics (the present 
plots were done with w\ = 10) and it is transformed at ground into a weight w, which is different 
for every muon, depending on its particular trajectory and energy. Weighted muons at ground were 
recorded and their distributions analyzed in comparison to the standard CORSIKA output. 



4.1. The energy distribution 

The energy at ground E f was histogramized and analyzed as a function of the impact point on 
ground (r, £). Typically, the muon energy is not directly measured by cosmic ray detectors since 
it would require extensive areas with particle detectors like those used in accelerator experiments 
and that is very expensive so far. Nevertheless, the spectrum of muons has an impact on other 
quantities that are measured by current air shower detector arrays, like the muon lateral distribution 
at ground, the arrival angle |26J], and the arrival time delay. 

Fig. Q2] displays a comparison between CORSIKA and the model for the normalized energy 
spectra of a deg and a 60 deg shower, at different distances from the shower core. The energy of 
muons decreases as ~ l/r and increases with the zenith angle 11141. Ilol. being the details given by 
the p t , z and E t distributions. Low energy muons dominate at large distances from the core. 
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Figure 12: Normalized energy spectrum of muons arriving at ground for a deg shower (left panel) and for a 60 deg 
shower (right panel), at different distances from the core as given by corsika compared to the prediction of the model. 



4.2. Apparent production depth distribution 

The production depth of the detected muons, the apparent MPD-distribution, follow distribu- 
tions whose shape changes with the observation coordinates. Muons can be understood as light 
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rays coming out of the shower axis with an radiance which is not isotropic, following an angular 
distribution given by the p t and Ej distributions as 



d 4 N 



dEdCldX 



1 cos a 

—h(X)f x (Ej, sin aEi)Ei— 

27r sin a 



(25) 



provided that the geomagnetic and multiple scattering effects are negligible. Since the angle a 
subtended from the observation point (r, changes with the production point at the shower axis z, 
the number of muons observed to come from a particular X will be different. Also, the production 
energy E t and the traveled distance I will be different and will induce different decay probabilities, 
modifying the shape of the total/true muon production depth distribution through p(E h I) into the 
apparent particular distribution of detected muons at that particular observation point. 



1 d 3 N 



We denote 



1 <fiN 
i- dXdrdl 



r dXdrd^ 
as dN/dX\(r t £) for short. 



-J 



d 4 N 



dEidQdX P 



p(E h l)dEi 



(26) 



• Corsika r=1 25 m 
r Corsika r=500 m 
a Corsika r=2000 m 

model 

Corsika h(X) 




X (g cm" 2 ) 



Figure 13: Comparison of several apparent MPD-distributions, dN/dX\( n 0, for a 40 deg shower at different distances 
from the core. The total/true MPD-distribution (h(X)) is also plotted for comparison. Normalizations are arbitrary. 



Fig. [T3]displays a comparison between CORSIKA and the present model of the apparent muon 
production depth distributions for a 40 deg shower at different distances from the core, where the 
distortions introduced in the dN/dX\^ distributions when compared to h(X) can be observed. 
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Fi gure 14: Comparison between model and CORSIKA of the normalized apparent MPD-distributions, dN/dX\^ r ^, 
for a 60 degree shower at 3 different distances from the core. The color histograms show the contribution of different 
energies. 



Fig. [14] displays a comparison at 60 degrees at different distances from the core. The distortions 
due to the observation point are small at high zenith angle. Different colors show the contributions 
from different energy of muons at ground. It can be seen that high energy muons tend to come 
from higher up in the atmosphere. At close distances to the core their contribution is enhanced. 

The dN/dXy,-^) distribution is never directly observed, but reconstructed from the arrival time 
or the arrival angle at ground. The correct inference of the total/true MPD-distribution, h(X), 
requires the knowledge of the exact dependence of dN/dX\ (r ^ with the observation point coordi- 
nates and detection energy threshold. dN/dX\^ explores different kinematic regions at produc- 
tion when reconstructed at different distances from the core. For instance, the algorithm proposed 



in fl 1 3h and 111 lh requires the conversion of each dN/dX\ {r ^ observed in each station to an universal 



distribution in order to sum up the contributions of all detectors in a single shower. 

4.3. Time distributions 

The total time delay is the sum of the different contributions calculated in the section 3. 

t = t g + t e + ts + tR em (27) 

where t g is the geometric delay, t £ is the kinematic delay, ts is the contribution produced by the 
geomagnetic field, and finally tR em accounts for the delay due to multiple scattering and the inacu- 
racies due to the approximations used. Fig. [15] displays the different delays for 3 different zenith 
angles, namely 0, 60 and 70 degrees, where the different contributions to the total delay are dis- 
played. The contribution to the multiple scatering is included in t Rem . At large distances from the 
core, the geometric delay is the most important. At distances typically from a few hundred meters 
to 1 km, the kinematic delay has a large impact. As we increase the zenith angle, the geometric 
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Figure 15: Average time delays for the different contributions to the total delay for 3 showers at 0, 60 and 70 degrees, 
from left to right. 




Figure 16: Comparison between the model and CORSIKA of the normalized time distributions for a 60 degree shower 
at 3 different distances from the core. The color histograms show the contribution of different energies. 



delay looses importance relative to total delay. At 500 m from the core, the geometric delay rep- 
resents - 70%,60% and 50% for 0, 60 and 70 degree showers, respectively. Fig. [Redisplays the 
overall time distributions at 3 distances of the shower core for a 60 deg shower. Filled histograms 
show the contributions of different muon energies at ground. High energy muons arrive earlier at 
ground. This is so because they are produced higher up in the atmosphere, and therefore have less 
geometric delay, but also because they have less kinematic delay. 
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The muon arrival time distributions can be used to extract important information. Far from 
the core, the time distributions are to a very good extent a one to one map of the apparent MPD- 
distributions. It can be determined by converting each muon time into a production distance, 
being the kinematic time a second order correction. Since the energy of each muon is typically 
not known, it is approximated by the mean value, taken from the energy spectrum at each obser- 
vation point as discussed in section |4"T1 and as it was explained in 11151 1 1411 . The energy would 
also determine the parameters of the multiple scattering delay distribution, although its concrete 
value follows a random distribution. The geomagnetic delay might take only two possible values 
depending on the charge of the muon. In general this technique will require a stringent r cut for 
those regions where the geometric delay is a large fraction of the total delay, in order to avoid 
distortions of the reconstructed dN/dX\^. An alternative method consists in fitting the time dis- 
tributions at once leaving a set of shape parameters on h(X) free. A detailed discussion of this 
procedure will be made elsewhere. 




Model E,>2 GeV 



Figure 17: Comparison between the model and CORSIKA of the muon lateral distribution at ground for 3 showers at 
0, 60 and 70 degrees, from left to right. The color histograms show the contribution of different energies. 



It is possible in principle to measure, or at least constrain, the shape of the muon spectrum. 
At distances close to the core the geometric delay is not dominant and the arrival time is mostly 
determined by the energy of each muon. A complete procedure to experimentally access this 
distribution is out of the scope of this paper, and will be analyzed elsewhere. 



4.4. Muon lateral distribution at ground 

The number of muons per surface area unit is p(r, £) = 



d 2 N 
rdrd( 



■ Fig. El displays examples of 



muon lateral distributions at ground for several zenith angles, where the contributions of different 
muon energies at ground were displayed in different colors. Low energy muons have a major 
impact on the fine details of the muon lateral distribution at ground. 
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Figure 18: Muon lateral distribution at ground for 3 different polar angles £ 



In vertical showers the number of muons per surface area does not depend much on £. As we 
increase the zenith angle, asymmetries appear because of the different propagation effects, mainly 
decay and geometry. The effects of the magnetic field become important above 60 degrees, and 
they completely dominate the distributions at very inclined showers, typically between 80 and 90 
degrees [24]. Fig. [T8]displays the muon density as a function of r for 3 different polar angles £ on 
a 70 deg shower. 

The shape of the ground distributions is fully determined by the distributions at production, 
h(X) and fx(Ej, p t ). A change in the overall muon content of the shower, No, produces a change in 
the muon density at ground, and therefore in the normalization of all distributions. The other main 
source of fluctuations comes from the depth of the first interaction, which directly affects h(X) by 
changing its maximum, X^ ax . The position of X^ ax directly influences all distribution at ground 
since it changes the total distance traveled by muons to ground. Fig. [19] left panel displays two 
different muon lateral distributions at ground for two different positions of X^ ax for a shower of 
60 degrees, where a change on the shape of the muon lateral distribution can be observed. Fig. [T9l 
right panel displays the normalized time distributions for (r, £)=(1000 m,-180 deg), for the same 
showers. 
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Figure 19: Muon lateral distribution (right panel) and normalized arrival time distribution at 2000 m and f = -180 
degrees (left panel) for two shower maxima X^ dx - 634 g cirT 2 , and X^ ax = 576 g cm" 2 



5. Average energy and transverse momentum distributions 

One of the main applications of the present model, is to be used in a global fit to extract in- 
formation on the total number of muons in the shower No, and the total/true production depth 
distribution, h(X), and its maximum, X^ ax . In order to do so, we must assume a fx(Ej, p t ) distribu- 
tion. 

The energy and transverse momentum distributions display more universal features when they 
are expresed in terms of X' = X - X^ ax , once the effects of the fluctuations induced by the first 
interaction point are removed. Fig. [31 right panel, displays an example of the energy spectrum for 
50 showers at X' = -300 g cm" 2 and X' = 300 g cm -2 . For comparison, Fig. @J left panel, displays 
the averaged distributions over the X' variable. The average energy and transverse momentum 
distributions do not change when changing the energy of the primaries (Figsj6]and|9l right panel), 
whereas they show mild differences between proton and iron primaries (Figs. [5]and|7l right panel), 
hadronic interaction models (Figs. [6]and[9l right panel). 

If we substitute fx(E t , cp t ) of a given shower by an average over showers of the same hadronic 
interaction model, primary, and zenith angle, < fx>(Ej, cp t ) >, and leaving only h(X) from the 
original shower, the ground density displays differences of about ~ 2% at 1000 m compared to 
the prediction it had occured if we used fx(Ei,cp t ), whereas the rest of the ground distributions 
remained unchanged. It is thus possible to use an universal energy and traverse momentum dis- 
tribution that depends only on X', where the position of X^ ax is naturally accounted for through 
X = X' + Xm ax . 

The systematics of any concrete aplication, including a global fit are to be studied and ac- 
counted for in each particular method and/or experimental setup. The effects of the choice of 
hadronic interaction model on < f x >(E h cp t ) > might introduce some systematics that should be 
also accounted for. 
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One could also think of a method to experimentally constraint the energy and transverse mo- 
mentum spectrum based on simultaneous observations of the ground distributions in different con- 
ditions. 

6. Conclusions 

In this paper we have built a model that explains the time delay, lateral distribution and depth 
profiles of muons at ground in terms of the distributions of energy, transverse momentum and 
production depth of all muons at the shower axis. The propagation of muons takes into account, 
continuous energy losses, decay, magnetic field effects and multiple scattering. The effects of 
multiple scattering in the time distributions, the lateral distribution, apparent depth profile, and the 
energy spectrum can be neglected at distances to the core larger than 100 m. The angular distribu- 
tion of muons at ground and the importance of the different effects will be analyzed elsewhere. 

This model can be used to experimentally reconstruct the distributions at production, in par- 
ticular the total/true muon production depth h(X) and the total number of muons No. The energy 
and transverse momentum distributions at production play a fundamental role on shaping the time 
distributions and muon lateral distribution at ground. The fx(E h p t ) distribution shows a quasi- 
universal behaviour when expressed in terms of X' = X-X^ ax and it can substituted by an average 
of showers performed over the X' variable in many applications. 

The model was implemented in a fast Monte Carlo, on two different steps: straight line prop- 
agation, and corrections, allowing the production of all distributions at ground starting from the 
distributions at production. 
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